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Histogram Monte-Carlo simulation results are presented 
for the magnetic-field - temperature phase diagram of the 
Ising model on a stacked triangular lattice with antiferromag- 
netic intraplane and ferromagnetic interplane interactions. 
Finite-size scaling results for this frustrated system at three 
points along the paramagnetic transition boundary are pre- 
sented which strongly suggest a line of triciritcal points at low 
field and a first-order transition line at higher fields. These 
results are compared with the corresponding phase diagrams 
from conventional mean-field theory as well as from the Monte 
Carlo mean-field calculations of Netz and Berker [Phys. Rev. 
Lett. 66, 377 (1991)]. 
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I. INTRODUCTION 

Despite being perhaps the simplest example of a ge- 
ometrically frustrated magnetic system, and having a 
long history of investigation,!"^ the Ising antiferromag- 
net on triangular and hexagonal lattioas continues to re- 
veal a number of remarkable fcatureai as well as elvdft 
concensus regarding some of its critical properties .ETlI 
The additional competition among interactions intro- 
duced by a magnetic field often gives rise to a wide va- 
riety of multicritical-point phenomena in frustrated sys- 
tems, with a concomitant change in critical behavior. Al- 
though mean-field theories usually can capture at least 
the qualitative structures of the magnetic phase diagrams 
as determined hv, more sophisticated methods, as well as 
by experiments,H there are notable exceptions, at least 
in the case of the XY modelJlj'EJ In this work, results 
of extensive Monte Carlo simulations of magnetic-field 
effects on the stacked triangular Ising model with an- 
tiferromagnetic intraplane and ferromagnetic interplane 
interactions are presented. This-, study was motivated 
by the work of Netz and Berkeiil3 (NB) who employed 
a "Monte-Carlo Mean-Field" theory to investigate the 
phase diagram. Although this method is essentially 
based on a mean-field approximation, it preserves the 
"hard-spin" condition that is believed to be essential 
for frustrated systems. In view of its relative simplic- 
ity, and the continuing controversy regarding the nature 



of phase transitions in frustrated antiferromagnets,! 
Monte-Carlo Mean-Field theory is a powerful technique 
which deserves further testing. 

The magnetic ordering process of the frustrated Ising 
antiferromagnet even at zero field has not been com- 
pletely resolved. As the temperature is lowered from the 
paramagnetic (P) phase, a transition to a partially dis- 
ordered period-three state (phase A) occurs where one- 
third of the sublattices remains paramagnetic. Modified 
mean-field theories predicts thaL,this phase remains sta- 
ble down to zero temperature.Ej A second peak in the 
specific heat at lower temperature observed in Monte 
Carlo (MC) results has been interpreted as being as- 
sociated with a transition to a ferrimagnetic state (phase 
B)B Such a transition also has been indicatedj-iji the 
results of NB. Other reports of MC simulations Ji3 how- 
ever, suggest that this peak is not associated with a phase 
transition but signals the onset of new types of excita- 
tions associated with the evolution of phase A into the 
zero temperature Wannier state.Q Notably, however, a 
low-temperature phase B is stabilized by the addition of 
next-nearest neighbor interactions, with an intermediate 
(between phases A and B) ferrimagnetic phase C, as de, 
termined by a variety of mean-field apprpyimat ionsr iM, 
and as recently observed experimentally.tJ 

Magnetic-field effects have been much less studied. 
The results of applying a modified mean-field theory to 
the triangular lattice is a phase diagram which exhibits 
a small region of phase A at low fields, a large region 
of phase B at higher field strengths, and with the Neel 
temperature being ja-multicritical point where phases A, 
B, and P coincide .Ii3 These results should also apply to 
the ferromagnetically coupled triangular layered system 
of interest here since mean-field theory is independent of 
space dimensions. The results of NB for the hexagonal 
system show two transitions at zero field, phase A sta- 
bilized in a small region at lower field values and higher 
temperatures, phase B occupying most of the phase di- 
agram, but with a possibly complex multicritical point 
occurring at a finite field along the paramagnetic phase 
boundary. An important aspect of these results is the 
observation that the P — B phase boundary is a first- 
order transition line due to the 3-state (Potts) symmetry 
of phase B. The transition at the Neel temperatureJa. 
believed to be of XY universality (due to frustration) BEI 



and thus, so should be the entire field-induced transition 
line A — P. NB also claim that the A — B transition 
is continuous. However, we have argued that these two 
states have the same symmetry so that the transition 
must be first order .t3 Other work on phase diagrams, of 
possibly indirect relevance, was based on a phenomeno- 
logical mpdel for the case of antiferromagnetic interlayer 
couplingJlj 

We study here the Ising Hamiltonian on a simple 
hexagonal lattice 



<ij> 



<kl> 



where J|| — —1 is the ferromagnetic interplane interac- 
tion, J± — 1 indicates the antiferromagnetic intraplane 
coupling which is frustrated for the triangular geometry, 
< i,j > and < k,l > represent near- neighbor sums along 
the hexagonal c axis and in the basal plane, respectively, 
with an applied field H . After a preliminary examination 
of the phase dijagram based on a simplified molecular- 
field calculation^ in the next section, the results of finite- 
size scaling of conventional and histogranxj MC data are 
presented. The utility of these approaches in determining 
phase diagraHiSiOf frustrated systems has been previously 
demonstrated.HO In contrast with other mean-field mod- 
els, our calculations yield a phase diagram where only 
phase B is present, even at zero field. Although the 
transition at TJ = is continuous, the P — B transition 
is found to be first-order for all H ^ {). The Neel tem- 
perature thus represents a tricritical point within this 
approximation. Interpretation of the MC data is more 
difficult due to the ambiguity in making a distinction be- 
tween phases A and B in the presence of a field, as well 
as the relatively large effects of critical fluctuations. It is 
clear, however, that phase A is the only state which oc- 
curs at _ff = (in agreement with other MC results and 
in contrast with NB). It is also evident that if this state 
exists at finite fields, it occupies only a very small region 
of the phase diagram, roughly consistent with the results 
oi NB. At high fields, the transition to the paramagnetic 
state is weakly first-order, in agreement with molecular- 
field theory and symmetry arguments of NB. At lower 
fields, however, the behavior appears close to tricritical 
for both of the points on the paramagnetic boundary that 
were examined. This is in contrast with the prediction of 
XY universality (NB) (as founda at H ^ 0) and suggests 
the possibility that phase A exists only at the paramag- 
netic boundary for H ^ 0. 



II. MOLECULAR-FIELD LANDAU THEORY 

Analyses of magnetic phase diagrams based on phe- 
nomenological Landau-type free energies for frustrated 
spin systems have proven to be quite successful in repro- 
ducing thcressential features of both MC and experimen- 
tal results. Usually, qualitative features can be captured 



using the more limited approach of deriving such a free 
energy from a molecular-field treatment of Hamiltonians 
such as that given by (1). This method allows for the 
understanding of relative phase stability from analytic 
expressions and makes symmetry arguments clear. For 
the present work, the free energy is expanded to sixth 
order in the spin density 



s{r) = 2[Xcos{Cl-r + (j))+Y] 



(2) 



where X is the amplitude of the modulated component, 
y is a uniform component, and Q = 47r/(3a)a is the 
period-three (in the basal plane) wave vector. Results for 
the free energy are similar to those of the phenomenolog- 
ical model given in Ref. |lj for the case of antiferromag- 
netic interlayer coupling, but with a few important dif- 
ferences. In the phenomenological model, the transition 
temperatures T^ and Ty, as well as the fourth and sixth- 
order coefficients, B and C, were treated as adjustable 
parameters. Within the molecular-field derivation, these 
quantities are determined by 

Tx = (-2J|| + 3Ji)/a, Ty = -(2J|| + 6Ji_)/a, (3) 

along with B — bT and C = cT, where the parameters a, 
6,|-apd c are given in terms of the total angular momentum 
jj23 For simplicity, we take here j = 1 so that, e.g., a — 3. 
In addition, there is a direct Zeeman coupling —HY, ab- 
sent in antiferromagnetic interlayer-coupling case of Ref. 

The three phases A, B, and C (labelled 1, 3, and 2, 
respectively in Ref. |l^) are conveniently characterised by 
Y and the phase angle 0: In phase A Y — 0, (f) — (2n + 
1)tt/6; in phase B Y ^ 0, (j) = mr/i; in phase C Y j^ 
0, and 4> takes on intermediate values. Stabilization of 
the 3-state Potts phase B is due to a term third-order 
in X, ~ Y X'^ cos{3(p) , and notably involves the unifopa 
component Y, even in the absence of an applied field.lij 

The resulting phase diagram is shown in Fig. 1. Only 
phase B occurs, even at zero field. The transition B — P 
is continuous at _ff = and is first order otherwise. The 
Neel transition thus represents a tricritical point. 

We note that within a corresponding phenomenologi- 
cal model, many different types of phase diagrams can be 
constructed, including those which feature all three types 
of ordered phases at zero field. For the purposes of com- 
parison with MC simulation results, only the molecular- 
field calculation is given here. 



III. PHASE DIAGRAM FROM MONTE CARLO 
SIMULATIONS 

Conventional Metropolis MC simulations were per- 
formed on the Hamiltonian (1) in an effort to determine 
the main features of the phase diagram. Periodic bound- 
ary conditions were used on LxLxL lattices with i=18- 
30. In some cases, a random initial spin configuration 



was used and thermodynamic averages were estimated 
after discarding the first 5 — 10 x 10^ MC steps for ther- 
malization. In most cases, however, the final well ther- 
malizcd configuration of a previous run was used for the 
initial spin directions. Averages were then made using 
5 — 10 X 10^ MC steps. The paramagnetic phase boundary 
was estimated by the location of extrema in a variety of 
thermodynamic functions, discussed in the next section. 
No such anomalies were observed at the A — B bound- 
ary line, which was crudely estimated at j&iite field by 
the calculated value of the phase angle (jyrl At H — 0, 
a well-defined signature of phase A is the absence of the 
uniform component Y. This quantity was found to ex- 
trapolate to zero at large L for all temperatures examined 
(T > 0.4) and no signature of a transition (T ^ 0.9) to a 
low-T phase B, as found in Ref. H, was observed. 

The phase angle cj) is not in general well defined since 
a multitude of values are energetically equivalent, e.g., 
(j) — 7r/6,7r/2, ... in phase A. MC results at iJ = in- 
dicated that an average over these values {{(j)) = 7r/4) 
was produced, suggesting a small energy barrier between 
these states. At finite field, however, a transition region 
(dependent on L and the number of MC steps used) 
from this averaged state to one with stabilized values of 
(f) — rnr/S was observed. We tentatively associate this 
observation with a possible A — B transition. Energy 
barriers between the degenerate states of phase B are 
likely higher due to the presence of the cubic Potts term 
discussed above and it thus appears that a particular 
phase angle is selected for a given MC run. Using MC 
data at L = 21, a region of possible stability for phase A 
determined in this way is indicated in Fig. 2. Note that 
the Y is always non-zero in the presence of a field due to 
the Zeeman interaction and thus is not a good indicator 
of an A — B transition. 

Although there is general agreement (except at H = 0) 
between the results of Fig. 2 and those of NB, the ex- 
istence of a finite region of the ordered state A in the 
phase diagram is suspect. We base this comment on the 
conspicuous absence of any anomalies in the thermody- 
namic functions, as well as the finite-size results for criti- 
cal exponents presented in the next section: These do not 
support the predicted XY universality of the boundary 
line where the A — P transition is supposed to occur. 

Additional structure found by NB in the region of the 
multicritical point where phases A, B, and P meet is 
beyond the scope of the present technique and thus could 
not be verified. 



IV. FINITE-SIZE SCALING 

The analysis of MC generated histograms used here to 
determine finite-siae scaling behavior is now well tested ia 
standard systemsE3 as well as in frustrated Heisenberg,E3 
XY,E3 and Ising models.Ll In particular, new results for 
the present model at H — Q yielded crtical-exponcnt 



estimates^ (3 = 0.341(4), 7 = 1.31(3), v = 0.662(9), giv- 
ing a — 0.012(30) by the Rushbrooke equality. These 
results agree well with estimates from standard calcula- 
tions on the unfrustrated Xy model, such as from renor- 
malization group analysis:EJ /3 — 0.349(4), 7 — 1.315(7), 
V = 0.671(5) (giving a = -0.013(15)). Finite-size seal, 
ing at first-order transitions is also-aow well understoodEj 
and tested in frustrated systemsJlj In addition, tricriti- 
cal behavior has been observed in exponent values for 
the frustrated XY modelO In an effort to determine the 
critical behavior of the paramagnetic phase boundary of 
Fig. 2, finite-size scaling at three points (H=0.5, 1.0, and 
3.0 as indicated on Fig. 2) was performed using the same 
quality of MC data as in Ref. ||. 

Scaling (of the form ~ L^) was performed at the crit- 
ical temperature, T^, estimated by the order-parameter 
cumulant (Um) crossing method^— as-|Used previously by 
us for frustrated spin systems B'EjE3 A variety of ther- 
modynamic functions were examined, including the spe- 
cific heat C, order parameter M, susceptibility x ^ 
(M2), energy cumulant U{T) = 1 - \{E^) / {E^ , and 
the logarithmic derivative of M, which is equivalent to 
Vi{T) = (ME)/{M) - {E), as well as the second loga- 
rithmic derivative of M, V2{T) = {M'^E)/{M^) - {E). 
MC runs were made on lattices i = 12 — 30, with aver- 
aging performed using 5 x 10^ MC steps for the smaller 
lattices, and up to 2 x 10^ MC steps for the larger lattices, 
after discarding the initial 1 — 5 x 10^ MC steps for ther- 
malization. Additional averaging was then made using 
the results from between 6 (smaller lattices) to 12 sepa- 
rate runs. This yields a reasonable 2.4 x 10^ MC steps 
used for averaging at the largest lattice size. Gross esti- 
mates of error bars, given as an indication of the general 
quality of the data, were made by taking the standard 
deviation of results from these runs. 



A. Scaling at H=0.5 

Results of the histogram analysis performed at iJ = 0.5 
and T = 2.99, (near the paramagnetic phase boundary) 
are presented here. Fig. 3 shows data at which Um{T) 
for L =12 and 15 cross corresponding results at larger 
values of L' , giving an estimate of the critical temper- 
ature T/v = 2.987(5). (This value compares with the 
result at zero field, T/v = 2. 9298(10). 0) Finite-size scal- 
ing of other thermodynamic functions was performed at 
T/v, as well as at T = 2.982 and T = 2.992, thus provid- 
ing a method with which to estimate errors associated 
with critical exponents.ES The extrapolated value of the 
energy cumulant, U* = 0.66664(3), shown in Fig. 4a, 
is close to the result 2/3 expected at a continuous tran- 
sition. Scaling of the other thermodynamic functions is 
shown in Figs. 4b-e. 

Critical exponent-ratios associated with these func- 
tions were estimated using all the data, as-well as with 
results at L = 12 and T = 15 excluded. CJ Except for 



the specific heat, only a small effect was observed in the 
estimated values. By this method, ajv was found to be 
0.85(4), 0.79(7), and 0.68(11), respectively, where quoted 
errors represent only the robustness of the fits. In con- 
trast, the same excercise yields 0.465(5), 0.455(4), and 
0.451(5) for (ijv. Results given on the figures are from 
fits made with only the L — VI data excluded. (Our 
conclusions are not affected by which set of data is used 
in the estimates.) The value Xjv — 1.947(5) was also 
determined from the data for V2 (not shown). Our fi- 
nal results at i/ = 0.5 are a = 0.40(13), /? = 0.23(5), 
7 = 1.09(12), and (taking the average of the two esti- 
mates for v) V — 0.51(3), where errors take account of 
uncertainties in Tjv. These values are not close to the ex- 
pected XY universality of an A — P transition (as found 
at if = 0). They are suggestive, however, of tricritical 
behavior, where a = 1/2, /3 = 1/4, 7 = 1, and v — 1/2, 



present caseJl^ In particular, finite-size effects appear to 
be especially important so that only results at the-Iarger 
lattice sizes may display the true critical behavior .EZl This 
is evident in the data shown for the energy cumulant, 
specific heat, and susceptibility (x) of Figs. 8a-c where 
volume-dependent scaling is observed only at larger L. 
Corresponding results for the susceptibility defined by 
x' ~ (M ^) — (Af )^ (considered to be less useful at a con- 
tinuous transition than x ^^ defined abova^), as well as 
V\ and V2, however, display volume dependence even at 
the smaller lattices, as seen in Figs. 8d-e. These behav- 
iors, along with the extrapolated value \J* — 0.66484(15) 
(considerably less than 2/3), are convincing indicators 
that the transition is indeed first order and involves the 
B phase. 



V. CONCLUSIONS 



B. Scaling at H=1.0 

Corresponding results at H — 1.0 (based on histograms 
generated at T = 3.07 for L = 12 - 27 and T = 3.068 for 
L = 30) are presented in Figs. 5 and 6. The critical tem- 
perature was estimated in this case to be Tjv = 3.067(2). 
Although the cumulant-crossing data for L=15 appears 
to be more widely scattered in Fig. 5 than in Fig. 3, note 
the difference in temperature scales by approximately 
a factor of two. The possibility that this transition is 
close to being first order is evidenced by the extrapo- 
lated result U* — 0.66656(8), which is slightly less than 
2/3 even within error. (As with all MC data, however, 
caution must be taken in assessing the reliabilty of error 
estimates.) Following the method outlined above, crit- 
ical exponents at this transition were estimated to be 
a = 0.59(14), /3 = 0.22(3), 7 = 0.97(8), and v = 0.46(2) 
(with the result l/i^ — 2.17(2) from V2 included). Again, 
these results are not those expected of XY universality, 
but (with the possible exception of i') are consistent with 
tricritical behavior. 



C. Scaling at H=3.0 

At the relatively high field value of ii — 3, a first- 
order transition from the paramagnetic to 3-state Potts 
phase B is expected. Significant scatter in the cumulant- 
crossing data of Fig. 7 is observed (from histograms gen- 
erated at T = 3.14 for L = 12 - 27 and T = 3.137 for 
L — 30), but again note the increased sensitivity of the 
temperature scale (a factor of three greater than that of 
Fig. 5). The transition temperature is estimated to be 
Tjv ~ 3.1365(10), with an error notably smaller than at 
the other two field values. 

Based on our previous experience with weak first-order 
transitions in frustrated spin systems, some care must be 
taken in the analysis and interpretation of the data in the 



Although general features of the phase diagram associ- 
ated with the ferromagnetically stacked triangular Ising 
antiferromagnet determined here by extensive MC sim- 
ulations are in agreement with the Mean-Field Monte 
Carlo results of Netz and Berker, some significant differ- 
ences are evident. In addition to finding only the par- 
tially ordered phase A stable at zero field, the results 
suggest that this phase is destabilized by any finite field. 
This suggestion is based primarily on finite-size scaling 
results at the supposed A — P transition where the ex- 
pected XY universality, found at if = 0, is not observed 
in the critical exponents. The results are, however, pos- 
sibly consistent with this boundary in the phase diagram 
being a line of tricritical points at which phases A, B 
and P are degenerate. We consider the possibility of non- 
universal critical behavior (as a function of field) unlikely 
as this phenomenon, as far as we know, has not previ- 
ously been observed in a three dimensional system. 

We find good evidence from finite-size scaling that the 
higher-field transition to the 3-state Potts phase B is in- 
deed first order. This result in fully consistent with N^B. 
and is also observed in the corresponding XY model.tiJ 
Somewhat surprising was the degree to which finite-size 
effects are important at this transition. For this reason, 
our original intention of studying this model with the 
addition of site disorder in an jeffort to observe a con- 
tinuous 3-state Potts transitionJ2j has been abandonded 
in favor of studying the same effect in the corresponding 
XY model.Eil It is hoped that the present work will serve 
to stimulate further theoretical effort as well as inspire a 
search for suitable experimentally testable realizations. 
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FIG. 1. Phase diagram from molecular- field theory. Re- 
gions labelled P and B represent paramagnetic and 3-state 
Potts phases, respectively. Broken line is a first-order phase 
boundary. 



FIG. 2. Possible phase diagram from Monte Carlo simu- 
lations, labelled as in Fig. 1, with the additional phase A 
representing the partially disordered state. Question mark 
indicates our uncertainty regarding the extent to which phase 
A exists in the phase diagram. Circles at H=0.5, 1.0, and 
3.0 represent field values at which finite-size scaling was per- 
formed. Squares denote possible locations of multicritical 
points. Broken and solid lines represent first-order and con- 
tirmous phase boundaries, respectively. 



FIG. 3. Results of applying the cumulant-crossing method 
to estimate the critical temperature at _ff = 0.5, where 
b=L'/L. 



FIG. 4. Finite-size scaling at the paramagnetic transition 
for H — 0.5 at Tjv = 2.987 where fitting is performed with the 
smallest lattice (L=12) omitted. Errors represent the robust- 
ness of the fit, except in the case of the energy cumulant (a), 
where it is determined by the uncertainty in Tjv. The specific 
heat, order parameter, susceptibility, and logarithmic deriva- 
tive of the order parameter are shown in (b)-(e), respectively. 



FIG. 5. Results of applying the cumulant-crossing method 
to estimate the critical temperature at _ff = 1.0. 



FIG. 6. Finite-size scaling at the paramagnetic transition 
for H = 1.0 at Tn = 3.067, as in Fig. 4. 



FIG. 7. Results of applying the cumulant-crossing method 
to estimate the critical temperature at H — 3.0. 



FIG. 8. Finite-size scaling at the paramagnetic transition 
for H — 3.0 at Tjv = 3.1365 with the assumption of volume 
dependence. 



